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Abstract 


V 

The  subject  of  this  thesis  was  an  attempt  to  find  a  periodic  solution  to 
the  equations  of  motion  of  a  high  eccentricity  (.7)  satellite  at  the  critical 
inclination,  with  period  equal  to  12  sidereal  hours,  and  with  the  earth's 
potential  represented  through  the  second  order  harmonics.  The  goals  of  the 
study  were  to  (1)  find  the  periodic  orbit  by  numerical  means,  (2)  examine  the 
stability  of  its  motion,  and  (3)  determine  the  characteristics  of  the  motion 
near  the  periodic  orbit  after  including  the  influences  of  the  sun  and  moon  in 
the  equations  of  motion. 

Success  was  obtained  in  the  two-body  case  for  all  orbits  and  in  the  two- 
body  +  J2  case  for  a  circular  equitorial  orbit  only.  In  the  two-body  +  J 2  case, 

no  periodic  orbit  could  be  found  with  non  zero  inclination  and/or  eccentricity 
t  0.  No  periodic  orbits  could  be  found  when  S22  and  C22  were  included  in  the 
equations  of  motion. 

^The  equations  of  motion  of  the  orbital  elements  were  examined  and  used 
to  explain  why  the  periodic  orbit  could  not  be  found. 


Determination  of  a  12-Hour  Periodic 
Orbit  at  the  Critical  Inclination 
with  the  Geopotential  Represented 
through  the  Second  Order 
Harmonics 

I  Introduction 

The  subject  of  his  thesis  is  an  attempt  to  find  a  periodic  solution 
to  the  equations  of  motion  of  a  high  eccentricity  (.7)  satellite  at  the 
critical  inclination  with  period  equal  to  12  sidereal  hours  and  with  the 
earth's  potential  represented  through  the  second  order  harmonics.  The 
orbit  is  shown  in  Figure  (1).  Apogee  is  in  the  northern  hemisphere  with 
argument  of  perigee,  u>,  equal  to  270°.  The  characteristics  just 
described  are  very  similar  to  those  of  the  Soviet  molynia  communications 
satellite  and  a  solution  would  be  directly  applicable  to  their  motion.  A 
solution  would  also  be  of  general  interest  because  of  the  special  charac¬ 
teristics  of  such  an  orbit.  The  12  sidereal  hour  period  and  the  critical 
inclination  both  correspond  to  special  cases  in  orbit  theory. 

The  equations  of  motion  of  a  satellite  about  the  earth  may  be  solved 
in  two  ways;  either  by  using  a  general  perturbations  method  to  obtain  an 
analytic  expression  for  the  motion  or  by  numerical  integration  to  obtain 
the  solution  over  some  particular  time  interval.  The  former  approach  is 
very  convenient  because  the  solution  can  be  used  for  any  orbit  after 
evaluating  the  initial  conditions  for  that  orbit.  However,  general 
perturbations  solutions  become  invalid  when  resonances  exist  between  two  o** 
more  frequencies  of  the  motion.  This  difficulty  is  known  as  "the  problem 
of  small  divisors"  (Ref  6:128-133). 

The  critical  inclination  and  12  sidereal  hour  orbits  are  special 
cases  because  they  correspond  to  resonance  conditions. 


Numerical  integration  of  the  equations  of  motion  is  very  straight¬ 
forward  and  does  not  suffer  from  resonance  problems.  However,  the  solution 
is  good  only  for  one  orbit  and  over  the  time  interval  through  which  it  is 
generated.  As  the  interval  of  the  integration  becomes  large,  the  amount 
of  labor  required  and  therefore  the  cost  of  generating  the  solution  become 
prohibative.  This  problem  can  be  solved  if  a  periodic  solution  can  be 
found.  It  has  to  be  generated  only  once  since  any  motion  outside  the 
solution  interval  is  simply  a  repeat  of  that  inside. 

Hori's  results  (Ref  7)  seemed  to  indicate  that  a  periodic  orbit  did  exist 
The  original  goals  of  this  study  became  (1)  find  the  periodic  orbit  by 
numerical  means,  (2)  examine  the  stability  of  its  motion,  and  (3)  determine 
the  characteristics  of  the  motion  near  the  periodic  orbit  after  including 
the  influences  of  the  sun  and  moon  in  the  equations  of  motion. 
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J~  II  Theory 
Equations  of  Motion 

The  equations  of  motion  of  the  satellite  are  determined  relative  to 
the  coordinate  system  shown  in  Figure  (2).  The  xyz-frame  is  fixed  to  the 
earth  which  is  rotating  with  respect  to  inertial  space  with  constant  angular 
velocity  flk.  <j>  is  the  angle  between  the  x-axis  and  the  projection  of  the 
satellite's  position  vector  in  the  x-y  (equitorial)  plane.  The  xyz-frame 
was  chosen  fixed  to  the  earth  so  that  later  the  Hamiltonian  of  the  motion 
would  be  constant.  The  e-frame  is  oriented  such  that  er  is  pointed  along 

A  A  A  A  A*  A 

the  outward  radial  ,  e.  =  kxe„,  and  e„  =  e.xe  .  This  set  of  coordinate 

<j)  r  o  <p  r 

A  A 

axes  becomes  undefined  when  er* k  =  +  1  so  that  the  equations  of  motion 
developed  using  these  coordinate  axes  cannot  be  used  for  polar  orbits. 

The  position  of  the  satellite  is 

R  =  rer  (1) 

where  r  is  the  distance  between  the  earth's  center  of  mass  and  the 
satellite  The  velocity  of  the  satellite  with  respect  to  inertial  space  is 

V  =  re„  +  ree.  +  rU+n)sinee,  (2) 

r  u  $ 

The  kinetic  energy  of  the  satellite  is 
T  =  JmV-V 

=  (r2  +  r2e2  +  r2($+n)Zsin2e)  (3) 

The  potential  energy  of  the  satellite  is  given  by  the  expression 
,  (Ref  2:275,  Ref  4:143) 
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V(r.e.*)  l  J  (V  p  "(cose) 

n=0  m=0 

[S„_sin  m  4>  +  Crt_cos  m  $]  (4) 

nm  nm 

where 

6  =  universal  gravitational  constant 
M  =  mass  of  earth 
Rg  =  radius  of  earth 

Pnm(cose)  =  associated  Legendre  polynomials 

Cnm’  snm  =  Geopotential  coefficients 

The  terms  with  m  =  0  are  called  zonal  harmonics  and  depend  on  latitude 

only.  The  coefficients  Cno  are  usually  written  Jn.  The  Sno  coefficients 

are  unimportant  since  sin *  0.  Terms  with  n  f  m  are  called  tesseral 

harmonics.  Sectorals  are  those  terms  with  n  =  m.  Tesserals  and  sectorals 

are  longitude  dependent  and  satellites  that  make  an  integral  number  of 

revolutions  about  the  earth  in  one  sidereal  day  are  in  resonance  with 

one  order  of  the  longitude  dependent  terms  of  the  geopotential.  A  twelve 

hour  satellite  is  in  resonance  with  the  n  *  2,  m  =  2  terms.  An  eight 

hour  satellite  is  in  resonance  with  the  n  >  3,  m  =  2  terms,  and  so  on. 

Because  the  12-hour  satellite  is  in  longidutinal  resonance,  the  C^2 , 

S2 2  coefficients  need  to  be  retained  in  the  potential. 

If  the  origin  is  chosen  to  coincide  with  the  earth's  center  of  mass 

terms  with  n  =  1  are  zero  (Ref  2:285},  the  J2  coefficient  is  the  most 

important  with  a  magnitude  on  the  order  of IxlO”3  compared  to  lxlO"6  or 

smaller  for  the  other  coefficients.  C2-|  and  S2-|  are  negligible  (Ref  8). 
o 

Pq  (cose)  *  1  so  if  CQ0  is  set  equal  to  1,  then  the  first  term  of  the 
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potential  is  recognized  as  the  two-body  potential.  Finally,  the  potential 
is  through  second  order 


V  = 


GMm 
r  ’ 


Gbin  R  2  , 

- ®  J,  4-(3cos2e  -  1) 

j«3  £  ^ 


.  Rg2  3sin2e  (S22sin2^  +  C22cos2$) 


(5) 


The  lagrangian  of  the  systems  is  given  by: 


L  *  T  -  V 


=  |-m  (r2  +  r2e2  +  r2(4>+ft)2sin2e)  +  ^ 


J2  Tj-  (3cos2e  -  1) 


ymRe2 

+ - — 

r3 


3sin2e(S22sin2i(>  +  C22cos29) 


(6) 


where  y(=GM)  is  the  gravitational  parameter.  The  Hamiltonian  of  the  systems 
is  given  by  H  =  tp^q^  -  L  where  q.  are  the  generalized  coordinates  and 
p^ (=3L/aq^ )  are  their  conjugate  momenta.  The  momenta  and  resulting  Hamiltonian 
are 


Pr  =  mr 
P6  =  mr2© 

pA  =  mr2(4>+n)$in2a 
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(7) 

(8) 
(9) 


\ 
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2  2 

H  *  tbi  [p  2  +  — —  +  — -  2mp$n]  -  ^jr 
2  r  r2  r2sin2e  r 

-  v”*e2_  J2  \  (3cos2e  -  1) 


ymR  2 

— —  3sin2e  (S22sin2<|>  +  C22cos24>) 


(10) 


Since  time,  t,  is  not  present  in  the  right  hand  side  of  Eq  (10),  the 
Hamiltonian  or  energy  of  the  system  is  constant.  For  two  body  motion  and 
motion  with  J2  not  equal  to  zero,  the  S22  =  C22  =  0  means  that  <f>  also  does 
not  appear  in  Eq  (10).  This  means  for  those  two  cases  p^  is  constant  along 
with  the  Hamilotnian.  The  equations  of  motion  are 


*i 


3H 

3Pj 


which  yield 


*1  =  *4 


y 


3 


*6 

y.j2sin2y2 


R 


y s  +  !§_  .1  -  1  j  1  (3cos2y2  -  l) 

4  v  Vsi'n2y2  V  V 


(11) 

(12) 

(13) 

(14) 


-  ^  sin2y2(S22sin2y3  +  C22cos2y3)  (15) 

y^ 
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*5 


y62c05y2  -30  1 
y-j2sin3yz  2  2  y^ 


s1n2y2 


+  ~  siny2cosy2  (S22sin2y3  +  C22cos2y3) 
yl 


(16) 


y6  =  —  s1"2y2  (S22COs2y3  -  C22s1n2y3} 
y1 


07) 


where 


yl  “r 


y2  a  e 


y,  *  ♦ 


y4“x 
y5  ’ir 

m 

*6 

m 


Units  of  mass,  time,  and  distance  were  chosen  such  that  Rg  =  1  and  u  =  1. 

Eqs  (12)  -  (17)  are  six  first  order  differential  equations  of  motion  for  the 
satellite. 


Resonant  Orbits 

All  two  body  orbits  are  periodic  with  respect  to  inertial  space.  In  the 
rotating  frame  only  resonant  orbits,  those  that  make  an  integral  number  of 
revolutions  about  the  earth  per  day,  are  periodic.  These  are  called  sub- 
sychronous  orbits  and  repeat  the  same  groundtrack  every  24  sidereal  hours. 
Although  they  have  different  inertial  periods,  each  has  the  same  period  in 
the  rotating  frame,  that  is,  24  hours.  Figure  (3)  shows  representative 
resonant  orbits  with  inertial  periods  of  8,  12,  and  24  sidereal  hours. 


Symmetry 

Figure  (3)  shows  that  for  the  12  hour  orbit,  apogee  occurs  at  two 
points  180°  apart.  At  apogee  r  *  0  so  that  the  orbit  intersects  the  x-axis 
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perpendicularly.  The  path  from  A'  to  A  (See  Figure  4)  must  be  Identical 
to  that  from  A  to  A*  since  they  represent  the  same  point  In  space.  Notice 
the  elliptical  shape  of  the  earth's  equator.  This  represents  the  presence 


of  the  longitude  dependent  terms  C22  and  S22  in  the  geopotential  and 
corresponds  to  the  location  of  the  Eurasian  and  American  continents  on  the 
globe.  By'choosing  apogee  to  line  up  along  either  the  major  or  minor  axis 
of  the  equitorial  ellipse,  the  paths  from  A  to  A'  and  vice  versa  become 
mirror  images  of  each  other.  The  orbit  is  symmetric  about  the  x-axis. 


Algorithm 

Consider  the  general  eauation  of  motion 


x  =  f(x(t),  t)  (18) 

Suppose  a  solution  xQ(t)  is  known.  How  does  a  nearby  solution  behave?  If 
at  the  time  tQ  the  solution  is  displaced  by  a  small  amount  dx(tQ),  then  the 
solution  at  a  later  time  can  be  expanded  as  a  Taylor  series  in  <sx(tg),  that 
is. 


if 


6xUq)  +  Higher  Order  Terms 


(19) 


6x(t)  ■  x(t)  -  xQ(t) 


(20) 


then 

«x(t)  =  $(t,t0)<5x(tg)  +  Higher  Order  Terms 


(21) 


where 


«  is  known  as  the  state  transition  matrix  and  has  the  property  «(tg,t0)  *  I. 
To  determing  how  «(t,tg)  propogates,  recall  from  Eq  (20) 


* 


X(t)  =  xQ(t)  +  6x(t)  (22) 

using  Eqs  (22)  and  (18) 

x(t)  =  f(xQ(t)  +  6x(t) ,t) 

=  f(x0(t),t)  +  |^y  6x(t)  +  Higher  Order  Terms  (23) 

xO 

or 

•  _ 

5x(t)  =  A  <sx(t)  +  Higher  Order  Terms  (24) 

where 


substituting  Eq  (21)  into  Eq  (24)  and  since  Sx(tQ)  is  arbitrary  gives 

•Ct.t0)  =  A(t)Kt,t0)  ’  (26) 

Eqs  (18)  and  (26)  can  be  integrated  at  the  same  time  giving  the  solution 
x(t)  and  its  local  variation  with  initial  conditions. 

Suppose  the  period  of  the  orbit  is  known  (or  chosen).  One  can  guess 
the  initial  conditions,  x0(0),  and  integrate  to  get  Xq(t),  Since  Xq(0) 
was  a  guess,  the  orbit  probably  will  not  close  and  some  closure  error  will 
exist. 

Ax  =  x(0)  -  x(t)  .  (27) 

Figure  (5)  shows  the  closure  error.  If  the  initial  conditions  are  changed 
by  an  amount  Sx(0),  then  the  solution  at  time  x  will  vary  by  an  amount 
6x(x).  The  orbit  will  still  not  necessarily  close  unless  ax  ■  6x(x)  -  6x(0). 
But 
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ix(x)  =  *(x,0)$x(0) 


1 

so  that  for  closure 

AX  *  *(x,0)6x(0)  -  6x(0) 

*  O(x,0)  -  l]sx(0)  (28) 

(See  Figure  (6))  Eq  (28)  represents  n  simultaneous  equations  in  the  n  unknowns 
6x(0) .  Since  higher  order  terms  were  neglected,  the  process  should  be 
iterated.  The  method  is  described  by  the  following  algorithm: 

1)  Guess  initial  conditions  Xq(0) 

2)  Integrate  to  get  xQCx)  and  $(x,0) 

3)  Solve  O(t,0)  -  I]6x(0)  *  Ax^,  where 

AXj  =  x^(0)  -  x^t) 

4)  x1+1(0)  =  ^(0)  +  6x(0) 

5)  Are  both  <5x(0)  and  ax^  small  enough? 

Yes  Stop 

No  -*•  Go  to  Step  2 

Singularity.  Solving  Eq  (28)  for  6x(0)  requires  the  calculation  of 
O(t,0)  -  I]"1.  For  a  periodic  Hamiltonian  system,  $(x,0)  will  have  2 
eigenvalues  equal  to  1  for  each  integral  of  the  motion  (Ref  9:430-433).  In 
the  two  body  and  two  body  +  J2  cases  there  are  two  integrals  of  the  motion, 
the  Hamiltonian  and  p^.  With  C22  and  S22  added,  p^  will  no  longer  be  constant 
but  the  Hamiltonian  will  remain  so.  This  means  that  [$(x,0)  -  I]  will  be 
singular  and  so  [*(x,0)  -  I]’1  is  indeterminate. 

To  avoid  the  singularity  problem,  hold  one  or  more  of  the  initial 
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I 


* 


1 


conditions  constant.  The  row  and  column  of  [ft(t,0)  -  relating  to  the 

constant  initial  condition  is  then  eliminated  resulting  in  a  matrix  whose 

dimension  is  smaller  than  *  and  is  not  singular. 

C22»  S22  motion  has  only  one  integral,  the  Hamiltonian.  The  magnitude 

of  the  2nd  order  sectoral  term  is  very  small  so  that  although  p.  is  no  longer 

ft 

constant,  the  change  in  p,  is  very  small.  This  can  result  in  the  same 

ft 

singularity  problem  described  above.  While  [ft(i,0)  -  I]  might  not  actually 
be  singular,  it  can  be  so  close  to  singular  that  the  computer  cannot 
distinguish  any  difference.  The  approach,  then,  is  the  same  as  before: 

Hold  one  of  the  initial  conditions  fixed. 

Modification.  By  choosing  the  initial  longitude  of  apogee  equal  to  0°, 
then  three  of  the  initial  conditions  become  known,  namely  r  =  e  *  ft  3  0. 

That  .leaves  three  unknown  initial  conditions  r,  e,  and  ft.  At  ft  =  180°, 
apogee  occurs  again  and  r  *  e  *  0.  The  modified  algorithm  proceeds  as  follows: 


1)  Guess  initial  conditions  r,  e,  ft 


f,  e,  and  ft  are  known 


2)  Integrate  to  get  Xq(t)  and  ft(r,Q) 

where  Xq(t)  and  ft(r,0)  are  defined  as  before 

3)  Solve  &'  (t.0)  -  I]sx  (0)  =  Ax.  (t) 

u  ki  ft 

where  x^  is  the  vector  of  known  initial  conditions  r  ,  xu  is  the 

e 

r 

unknown  initial  conditions  e  ,  and$'  (x,0)  is  the  submatrix  of 

ft 

ftCt.O)  defined  by  *'(t,0)  =  3xk(T> 

axJPT 


4)  x  (0)  =  x  (0)  +  5x  (0) 
u1+l  ui  u 
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5)  Evaluate  the  residuals  ax.,  and  6xu(0) 
Are  both  small  enough? 

Yes  Stop 

No  Go  to  Step  2 


Program 

The  integration  of  Eqs  (12)  -  (17)  and  (26)  was  carried  out  using  a 
numerical  integration  package  named  ODE  (Ordinary  Differential  Equations). 
This  package  is  part  of  the  CC6600  library  for  use  on  the  CDC6600  computers 
at  AFIT.  The  program  is  listed  in  Appendix  B. 

The  elements  of  the  matrix  A  were  needed  in  the  program. 


(29) 


The  partial  derivatives  of  Eqs  (12)  -  (17)  were  taken  and  elements  of  A 
determined  by  Eq  (29).  Appendix  A  lists  the  derivatives. 


Verification  of  A  and  $  Matrices 
Recall  Eqs  (25),  (26),  and  (29) 

A  =  — 

**  (25) 

Ht.tg)  -  A(t)  $(t,tQ)  (26) 


aij  •  (29 

,xj 

and  that  $(0,0)  *  I.  In  order  to  verity  that  the  partial s  3f./ax.  were 

*  3 

taken  correctly  and  that  $  is  being  propogated  correctly,  alternate  methods 
were  used  to  form  both  A  and  $  and  compared  to  the  results  given  by  Eqs 
(25)  and  (26). 
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Elements  of  A.  Given  a  solution  x(t),  then  from  Eqs  (12)  -  (17)  f(t) 
can  be  evaluated.  Perturbing  the  jth  element  of  x(t)  by  an  amount  ax.  gives 

J 

the  solution  xp(t).  fp(t)  can  then  be  evaluated.  The  jth  column  of  A  is  then 
approximately 


(30) 


Each  column  of  A  can  be  obtained  in  a  similar  manner.  The  results  of  Eq  (30) 
can  be  compared  to  those  of  Eq  (25)  and  checked  for  agreement. 

$  Propagation.  The  proof  that  *  is  being  propogated  correctly  is  very 
similar  to  the  one  just  discussed  for  the  A  matrix.  Given  the  initial  conditions 
xQ(0)  integrate  to  get  x(t).  Next  perturb  the  jth  element  of  xQ(0)  to  get 
xp(0).  Integrate  again  to  get  xp(t).  The  jth  column  of  *  is  approximately 


♦.  (t.0)  =  "xP(t)  “ 

Pj  AXj (0) 


(31) 


Each  column  of  $p  can  be  generated  in  a  similar  manner.  The  results  of 

Eq  (31)  and  the  integration  of  Eq  (22)  can  be  compared  and  checked  for  agreement. 
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Ill  Results  and  Discussion 


Results 

The  formation  of  the  A  matrix  and  propogation  of  the  $  matrix  were 
verified  using  the  methods  described  earlier.  Tables  (1)  and  (2)  list  the 
results  of  the  comparisons  showing  good  agreement  in  both  cases.  Very  large 
values  were  used  for  J2,  S22>  and  C22  50  arW  errors  in  the  terms 
associated  with  those  coefficients  would  be  discernable. 

Another  check  of  the  program  consisted  of  using  the  equations  of  motion 
(Eqs  (12)  -  (17))  to  verify  the  location  of  four  equilibrium  points  aligned 
along  the  principal  axes  of  the  equitorial  eclipse  (Ref  5:167-169).  Setting 
the  right  hand  sides  of  Eqs  (12)  -  (17)  to  zero  and  with  e  =  90°  gave  the 
four  equilibrium  points. 


With  J2  *  S22  =  C22  *  0  the  program  was  tested  to  see  if  two  body 
resonant  orbits  converged.  First,  precalculated  two  body  elements  were 
inputted  into  the  program  and  tested  for  immediate  convergence.  Convergence 
was  obtained  in  one  iteration  in  most  instances  and  always  within  two  itera¬ 


tions.  Second,  .elements  that  were  slightly  displaced  from  the  two  body  elements 
were  inputted  into  the  program.  Convergence  to  the  two  body  orbit  was  obtained 
very  quickly.  This  was  done  for  both  circular  and  elliptical  orbits,  inclined 
to  and  in  the  plane  of  the  equator.  Last,  to  show  that  non-resonant  orbits 
do  not  close  in  the  rotating  frame,  an  orbit  was  inputted  into  the  program 
which  then  tried  to  force  convergence  to  the  non-resonant  orbit.  Convergence 
could  not  be  obtained  which  was  interpreted  as  proof  that  non-resonant  orbits 
are  not  periodic  in  the  rotating  frame. 

It  can  be  easily  seen  that  the  modified  algorithm  has  many  variations 


depending  on  the  choice  of  known  and  unknown  initial  conditions.  For  the 
two  body  problem,  if  r^(0)  and  t  are  known,  then  $(0)  can  be  calculated. 
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TA9L7  i.  RESULTS  Or  r  MATRIX  V7R1 rI7ATT0N 
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TA3L7  2.  RESULTS  OF  AHI  MATRIX  VERIFICATION 


PHI  MATRIX 
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Check  of  a  matrix  had  R  =  e  =  <f>  =  P r  =  Pe  =  P<j>  =  1 ,  and  time  equal  zero. 
Check  of  PHI  matrix  had  initial  conditions  R  -  7.0796,  e  =  26.6°,  $  =  22.5° 
P  =  Pe  =  0,  Pa  =  .6525,  and  was  integrated  to  t  =  5 3,3979.  Both  programs 
and  AXj  =  l.E-5  and  Cg2  =  $22  =  ^2  =  ^*083E-2. 


•  •  ay  PKACTI^BLB 


19 


Since  the  integration  begins  at  apogee  r(0)  and  e(0)  are  known.  e(0)  and 

<j>(0)  are  arbitrary  and  so  are  also  known.  The  algorithm  is  modified  so  that 

r 

e 

xk(x)  =  $  .  xutO)  =  {p^} 

r 

e 

next,  J2  was  included  in  the  equations  of  motion  and  the  same  choice  of  known 
and  unknown  initial  conditions  was  made.  Convergence  was  obtained  for 
circular,  equitorial  orbits  but  periodic  orbits  could  not  be  obtained  that 
were  inclined  to  the  equator  or  elliptical.  When  C22  and  S22  were  included 
along  with  J2  then  convergence  could  not  be  obtained  in  any  case. 

Discussion 

After  repeated  failures  to  find  the  periodic  orbit,  it  became  obvious 
that  the  algorithm  being  used  was  not  working.  Either  the  computer  algorithm 
and/or  equations  of  motion  were  incorrect  or  the  periodic  orbit  being  searched 
for  did  not  exist.  At  first  the  former  was  thought  to  be  the  most  likely, 
but  the  good  agreement  obtained  in  the  check  of  the  A  and  $  matrices  along 
with  the  success  of  the  two  body  orbit  check  indicated  that  the  program  was 
working  correctly.  If  the  periodic  orbit  does  not  exist,  then  the  question 
becomes  "why  not?" 

This  study  was  begun  on  the  basis  of  Hori's  (Ref  7)  results  which  were 
believed  to  indicate  that  the  periodic  orbit  did  exist.  It  is  now  believed 
that  those  results  were  misunderstood.  General  perturbation  theory  gives 
for  the  time  derivatives  of  the  orbital  elements  (Ref  3) 


0-5  cos2i ) 


(33) 


i  =  Inclination 
a  =  semi-major  axis 
n  =  longitude  of  the  ascending  node 
u  =  argument  of  perigee 
e  =  eccentricity 
n  =  mean  motion 

At  the  critical  inclination  u  *  0  (the  critical  inclination  is  known  as  such 
because  it  is  the  inclination  at  which  the  rotation  of  the  apses  changes 
direction).  Eq  (32)  predicts  the  rotation  of  the  node  and  could  explain  why 
the  algorithm  failed  to  find  a  periodic  orbit.  At  first  it  was  thought  that 
since  the  node  seemed  to  rotate  by  180°  each  period  that  the  very  small  rotation 
predicted  by  Eq  (32)  would  be  negligible.  This  may  not  be  the  case.  The 
success  obtained  in  the  two  body  +  Jg  case  for  the  circular  equitorial  orbit 
could  have  disguised  the  nodal  regression  which  cannot  be  seen  in  such  an 
orbit.  If  nodal  regression  provides  the  clue  as  to  failure  occurred,  then 
the  supposed  period  of  12  hours  may  be  incorrect.  The  correct  period  of  the 
periodic  orbit  if  it  exists  is  probably  related  to  the  magnitude  of  Jg.  In 
fact,  Barrar's  work  (Ref  1)  proves  the  existence  of  periodic  orbits  at  the 
critical  inclination  about  an  oblate  earth  whose  periods  are  on  the  order  of 

1  1 
o23/2 
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Appendix  A 
Elements  of  A  Matrix 


*11  *  *12  ‘  *13  '  *15  "  *16  “  0 


*14  *  1 


*21  *  -2  X 


a22  *  a23  =  a24 s  a26  *  0 


a  =  "2y6 
a31  — - 


CA-1) 


CA-2) 


CA-3) 


CA-4) 


(A-5) 


(A-6) 


y12sin3yz 


_  -2yfi  cosy- 

a32  '  — ^ - 2_ 

y^sin3^ 


(A-7) 


a33  =  a34  =  a35  =  0 


(A-8) 


y^z  sin  2y2 


(A-9) 


.  =  -3y  2  3yfi2  .  9  .  6d9 

41  y-if  y-i r  y-3  (3cOS2y2  -  1) 

yl  yl  sin2y2  yl  yl  Z 


+  sin2y2  (S22  sin  2y3  +  C22  cos  2y3) 


(A-10) 
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*42  *  -y  c°sy2  +  a> 

y13sin3y2  y1 


cosy2  sin  y2 


—  siny2  cosy 2  (S22sin2y3  +  C22cos2y3) 

j'i4 


a43  "  —  Sitl2^  (S,-COS2y,  -  CjoSinSyg) 


a44  =  0 


a45  *  — ^ 


*1 


a46  ’  2/6 


y^sin2 


a51  = "  2y62cosy2  +  9  ^2  sin2y2 
y^3sin3y2  2  y-,4 


-  18 

^  siny2cosy2  (S22sin2y3  +  C22cos2y3) 


Hz  =  -  3Vcos2*z  -  V 


3J, 


y^sin^g  y^sinyg  y1 


-  cos2y2 


+  -  (cos2y2  -  sin2y2)  (S22sin2y3  +  C22cos2y3) 
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a53  ~  —  sir\y2cosy2  (S22cos2y3  -  C22sin2y3) 


(A-11 ) 

(A-12) 

(A-13) 

(A-14) 

(A-15) 

( A- 16) 


(A-T7) 

(A-18) 
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*54  ’  *55  ‘  0 

(A-19) 

.«  -  2y6C0Sy2 
y^sin^g 

(A-20) 

a61  =  sinZy2  (S22cos^y3  "  C22sin2y3) 

1 

(A-21 ) 

ag2  =  —  siny2cosy2  (S22cos2y3  -  C22sin2y3) 
yl3 

(A-22) 

ag3  =  "  —  sinzy2  (S22sin2y3  +  C22cos2y3) 

Y1 3 

(A-23) 

3 64  =  a65  =  a66  =  0 

(A-24) 
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APPENDIX  B 

Computer  Program  to  Calculate  Periodic  Orbits 
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